Computer aided treatment planning and visualization with image registration and fusion

ABSTRACT

A computer based system and method of visualizing a region using multiple image data sets is provided. The method includes acquiring first volumetric image data of a region and acquiring at least second volumetric image data of the region. The first image data is generally selected such that the structural features of the region are readily visualized. At least one control point is determined in the region using an identifiable structural characteristic discernable in the first volumetric image data. The at least one control point is also located in the at least second image data of the region such that the first image data and the at least second image data can be registered to one another using the at least one control point. Once the image data sets are registered, the registered first image data and at least second image data can be fused into a common display data set. The multiple image data sets have different and complimentary information to differentiate the structures and the functions in the region such that image segmentation algorithms and user interactive editing tools can be applied to obtain 3d spatial relations of the components in the region. Methods to correct spatial inhomogeneity in MR image data is also provided.

FIELD OF THE INVENTION

[0001] The present invention relates generally to surgical planning andmore particularly relates to systems and methods of registering,segmenting and analyzing image data for three dimensional interactivecomputer visualization and surgical planning and follow-up evaluation.

BACKGROUND OF THE INVENTION

[0002] In many areas of medical treatment, it would be beneficial for amedical practitioner to be able to visualize a region for whichtreatment is contemplated and to accurately simulate the contemplatedtreatment. By visualizing the effect of the simulated treatment andaltering the proposed treatment to optimize the results in a virtualsetting, results can be improved and risks associated with the actualtreatment can be reduced. This is particularly true in the case ofinvasive procedures such as surgery, biopsies and prosthesisimplantation. The virtual setting would serve both as a tool for theguidance for actual treatment and as a “gold standard” for evaluation ofthe actual treatment and for follow-up management. Such a system canalso provide an intuitive tool for medical training.

[0003] Preoperative imaging, such as by computerized tomography (CT),magnetic resonance imaging (MRI), ultrasound (US) and the like, isconsidered an important element in surgical planning. These conventionalimaging technologies provide two dimensional (2D) image slices whichillustrate key anatomic structures in a region of interest. However, the2D slice images are limited in their ability to represent the spatialrelationships between adjacent and interrelated structures. For example,a retro-displaced bifurcation point of carotid arteries might makestenosis analysis and plaque removal planning difficult.Three-dimensional (3D) information would be very helpful for planningsurgical therapy.

[0004] Traditionally, spatial anatomical relationships could only besurmised by mentally integrating sequential 2D slice images. However,the advent of 3D computer graphical visualization techniques andhigh-resolution image scanning now allow 3D images (either surface-basedor volume-based) to be constructed from sequential slice images and bedisplayed on a computer screen. Three-dimensional relationships betweenadjacent organs can be shown by interactively manipulating these virtualorgans on the screen using a mouse or some other interactive devices.Over the past decade, many applications of such techniques in a numberof areas of medicine, including otology, and neurology have beenexplored.

[0005] To achieve more accurate tissue classification and insightregarding tissue functionality, a series of single modality images andmulti-modality images can be utilized. The single modality images can beof a common region acquired at different times or differentorientations. In multi-modality images, the same region is imaged usingtwo or more imaging technologies, such as CT, MRI and US in order totake advantage of the properties of each technology since a certainmodality image may be sensitive to certain kind of tissues or theirfunctions. By using single modality image series and multi-modalityimages, more information can be obtained, and limitations of certainmodality can be mitigated.

[0006] The use of multiple image sets (including multi-modality imagesand or single modality image series) would be useful to perform virtualvisualization and treatment planning for carotid artery stenosis andother conditions. Carotid artery stenosis is the most common cause ofstroke, which is a major health care problem, that affects more than700,000 Americans each year. In fact, this condition is the leadingcause of disability and the third leading cause of deaths in the UnitedStates, after cardiovascular disease and cancer. Stenosis arises fromthe formation of plaque (consisting of calcification, cholesterol, andblood elements). The surgical procedure to remove plaque from a neckartery is called carotid endarterectomy. Prior to endarterectomy, thedegree of stenosis needs to be measured and the position of the plaquemust be localized. Currently available methods for evaluating carotidstenosis include, for example, carotid Doppler ultrasound and contrastangiography, which have been demonstrated for accurate determination ofthe degree of luminal stenosis. However, luminal narrowing is anindirect marker of plaque size and may underestimate the plaque burdenas plaque may grow inside the vessel wall from the lumen towards theouter surrounding tissues. It is desirable that plaque size (bothinwards and outwards boundary) and composition are accurately measured.These measures are of value since both plaque rupture and plaque removalcarry risk. The measures relating to plaque composition and thelikelihood of rupture can offer valuable risk assessment for thedecision whether or not to proceed with a plaque removal operation.

[0007] Currently, MR data, CT data and US data each provide some insightinto the structure, nature and position of stenosis within an artery.However, neither imaging technology alone is sufficient to provide acomplete analysis of the size, composition and position of the plaquebuildup. The benefits of these imaging technologies are largelycomplimentary and there would be a great benefit in having the abilityto readily register images of a region using these complimentarytechnologies and to fuse the images into a single display. There is alsoa benefit in having the ability to register and fuse a series of singlemode images to capture complimentary information contained therein.

[0008] In MR image data, the inter-slice and even intra slice image dataoften contains spatial inhomogeneity which adversely affects automatedanalysis of these images. Correction for spatial intensityinhomogeneities in inter- and intra-slices provides improved results inautomated quantitative analysis of MR images. The inhomogeneitieslocally alter image intensity mean, median, and variance. It is knownthat they vary over time and with different acquisition parameters andsubjects. For virtual treatment planning, a general purpose correctionalgorithm that can work for any scanning protocol and anatomy region isdesirable. The concept of correcting image artifacts can be applied tothe beam-hardening problem encountered in CT images as well as theattenuation distortion in US images.

[0009] Accordingly, there remains a need for improved medical treatmentplanning tools for optimizing the procedures and for evaluating theresults.

SUMMARY OF THE INVENTION

[0010] It is an object of the present invention to provide a method offlexible image registration suitable for use in three dimensional (3D)visualization systems.

[0011] It is a further object of the present invention to provide amethod of visualizing a region using multiple image sets, such as a setof single-modality images or multi-modality images, in a manner wherecomplimentary image data properties can be fused in a common, registereddisplay.

[0012] In accordance with one embodiment of the present method, acomputer based method of visualizing a region using multiple image datasets is provided. The method includes acquiring first volumetric imagedata of a region and acquiring at least second volumetric image data ofthe region. The first image data is generally selected such that thestructural features of the region are readily visualized. At least onecontrol point is determined in the region using an identifiablestructural characteristic discernable in the first volumetric imagedata. The at least one control point is also located in the at leastsecond image data of the region such that the first image data and theat least second image data can be registered to one another using the atleast one control point. Once the image data sets are registered, theregistered first image data and at least second image data can be fusedinto a common display data set.

[0013] The method can be used in connection with the visualization ofbranched structures, such as arterial structures. In this case, abifurcation point in the arterial structure can be selected as thecontrol point.

[0014] The first image data and the additional sets of image data can beof the same region at different times, of the same region with differentpatient orientations, or of the same region using different imagingmodes. The various sets of images, either single-modality image sets ormulti-modality images, are selected to provide complimentary imagingcharacteristics. For example, in the case of visualizing an arterialstructure for carotid stenosis, a first imaging mode can be used todelineate the arterial lumen, a second imaging mode can be used todelineate fatty components of plaque deposits and third imaging mode canbe used to delineate calcification components of plaque deposits.

[0015] A method is also provided for correcting spatial inhomogenietieswhich may occur in MR image data. After acquiring the MR image data, abias field associated with the image data is estimated. The bias fieldis then applied to the MRI data to correct for the spatialinhomogenieties.

[0016] A method is also provided to assist a user, such as a physician,in editing the computer-aided, automatically processed results of avisualization, such as the segmentation of the tissues within the regionfrom the multiple image sets. The user employs a graphical userinterface to introduce editing to the multiple image sets. Such editingcan also be used to manually optimize the registration and fusionoperations for the multiple image sets.

BRIEF DESCRIPTION OF THE DRAWING

[0017] Further objects, features and advantages of the invention willbecome apparent from the following detailed description taken inconjunction with the accompanying figures showing illustrativeembodiments of the invention, in which:

[0018]FIG. 1 is a simplified flow diagram illustrating an overview of amethod for computer aided treatment planning and interactivevisualization;

[0019]FIG. 2 is a simplified block diagram of a system suitable forperforming the present method of computer aided treatment planning andinteractive visualization;

[0020]FIG. 3 is a flow diagram illustrating an overview of a flexibleregistration process for aligning multiple sets of image data, such asin multi-modality imaging and single-modality image series acquisition;

[0021]FIG. 4 is a pictorial diagram illustrating a section of a bloodvessel with plaque build up resulting in stenosis;

[0022]FIG. 5 is a pictorial diagram illustrating a branched section of ablood vessel with plaque build up resulting in stenosis, the branchedsection featuring a bifurcation point which can be used as a controlpoint for image data registration;

[0023]FIG. 6 is a pictorial diagram of a branched section of a bloodvessel further illustrating the location of the bifurcation point whichcan be used as a control point for image data registration;

[0024]FIG. 7A further illustrates the determination of the control pointfor a spherical model;

[0025]FIG. 7B illustrates the application of the control point andcoordinate system of FIG. 7A as applied to a deformable model.

[0026]FIG. 8 is a flow diagram further illustrating the method of FIG. 3as applied to a method of visualizing and analyzing carotid stenosis.

[0027]FIG. 9 is a flow diagram illustrating a process for correctingspatial inhomogenieties in MRI image data.

[0028]FIG. 10 is a pictorial diagram illustrating an exemplary graphicaluser interface display illustrating segmented multi-modality imagesafter a fusion operation.

[0029] Throughout the figures, the same reference numerals andcharacters, unless otherwise stated, are used to denote like features,elements, components or portions of the illustrated embodiments.Moreover, while the subject invention will now be described in detailwith reference to the figures, it is done so in connection with theillustrative embodiments. It is intended that changes and modificationscan be made to the described embodiments without departing from the truescope and spirit of the subject invention as defined by the appendedclaims.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0030]FIG. 1 is a flow chart which illustrates an overview of thepresent method of computer aided treatment planning and interactivevisualization which is generally performed on a computer based system,such as that illustrated in FIG. 2. The invention will be described interms of medical applications performed on human patients and in thecontext of medical treatment, such as examination ( e.g., measuringcarotid stenosis and quantifying plaque components), surgery (e.g.,removal of carotid plaque), prosthesis implantation, biopsy, medication,therapeutic radiation, therapeutic ultrasound and the like. It will beappreciated, however, that the invention is not limited to humanpatients, nor to the exemplary list of treatments referenced. The termtreatment is used to mean an intervention in a region, such as but notlimited to tissue, that is intended to effect an alteration of theregion.

[0031] Referring to FIG. 1, the method includes the initial step ofgenerating a three dimensional (3D) image representation of a region forwhich some form of medical treatment or intervention is contemplated(step 102). Generating such a 3D image representation generally involvesacquiring a sequential series of 2D slice images, such as from a spiralcomputed tomography (CT) scanner, magnetic resonance imaging (MRI)scanner or ultrasound scanner (US) and transforming this 2D image datainto a volumetric data set which provides a 3D representation of theregion on a 2D display, such as a computer monitor. Such a technique iswell known in the art, and is discussed, for example in U.S. Pat. No.5,971,767 to Kaufman et. al., which is hereby incorporated by referencein its entirety.

[0032] After the 3D image is presented to a user, such as a physician,some form of virtual intervention, which simulates at least a portion ofa proposed treatment, is applied to the 3D image (step 104). The virtualintervention can take on several forms, such as the quantifying ofarterial stenosis and plaque components, removal of tissue or arteryplaques, the repair or reconstruction of a diseased or malformed organ(e.g., inner ear, lungs, liver, joints, etc.), the placement of aprosthetic implant, the placement of a stent graft, the placement ofbiopsy needle, the placement of therapeutic radiation and the like.

[0033] Using the resulting 3D image, and possibly the assistance ofcomputer generated models of the applied intervention, the results ofthe virtual intervention can be evaluated and warnings can be generatedindicative of high levels of risk attendant with the proposedintervention (step 106). Based on the displayed results, and anywarnings provided, the user can repeatedly modify the proposedintervention (step 108) and apply the modified intervention to the 3Dimage (step 104) until a satisfactory result is ascertained or it isdetermined that the proposed treatment is not feasible. Severalalternative interventions can be saved in a database to compare therisks and efficacy of proposed alternative intervention plans. Thequantified arterial stenosis and plaque components will provide ameasure of risk factor for plaque rupture, which can result in stroke,pulmonary embolism, and occlusions in the liver and kidneys.

[0034] After the proposed intervention is finalized, the finalintervention can be simulated and the results fully applied to the 3Dimage (step 110). The user can then view the results and navigate in andaround the region to determine the efficacy of the proposed treatment(step 112). The planned results can then be used as a guide for theactual treatment with coordinate registration between the virtual modeland the patient and as a gold standard to evaluate the actualintervention during post-intervention follow up examinations.

[0035]FIG. 2 is a simplified diagram of an exemplary system forperforming the present computer aided treatment planning and interactivevisualization methods. In this exemplary embodiment, a patient 201 liesdown on a platform 202 while scanning device 205 scans the area thatcontains the organ or organs which are to be examined. The scanningdevice 205 contains a scanning portion 203 which acquires image data ofthe patient and an electronics portion 206. Electronics portion 206generally includes an interface 207, a central processing unit 209, amemory 211 for temporarily storing the scanning data, and a secondinterface 213 for sending data to the virtual navigation platform.Interface 207 and 213 can be included in a single interface component orcan even be the same component. The various operational components andsubsystems in electronics portion 206 are connected together withconventional connectors.

[0036] The data from the scanning portion 203 is generally in the formof a stack of two dimensional image slices of a region of interest,which are provided from conventional spiral CT, MRI or US scanners.Central processing unit 209 converts the scanned 2D data to a 3D voxeldata representation, in a manner known in the art, and stores theresults in another portion of memory 211. Alternatively, the converteddata can also be directly sent to interface unit 213 to be transferredto the virtual navigation terminal 216. The conversion of the 2D datacould also take place at the virtual navigation terminal 216 after beingtransmitted from interface 213. Preferably, the converted data istransmitted over carrier 214 to the virtual navigation terminal 216 inorder for an operator to perform the computer aided treatment planning.The data can also be transported in other conventional ways such asstoring the data on a storage medium and physically transporting it toterminal 216 or by using satellite transmissions.

[0037] The scanned data need not be converted to its 3D representationuntil the visualization rendering engine requires it to be in 3D form.This may save computational steps and memory storage space.

[0038] Virtual interactive visualization terminal 216 includes a screen217 for viewing the image data, an electronics portion 215 and interfacedevice 219 such as a keyboard, mouse or track ball. The electronicsportion 215 generally includes a interface port 221, a centralprocessing unit 223, other components 227 necessary to run the terminaland a memory 225. The components in terminal 216 are connected togetherwith conventional connectors. The converted voxel data is received ininterface port 221 and stored in memory 225. The central processor unit223 then assembles the 3D voxels into a virtual representation which canbe displayed on screen 217. Preferably, a graphics accelerator which isoptimized for volume rendering can also be used in generating therepresentations. The virtual interactive visualization terminal 216 canbe embodied using a high speed graphics work station, such asmanufactured by Silicon Graphics, Inc., or in a high speed personalcomputer, such as an IBM compatible computer with a Pentium III (orhigher) processor having a 1 GHZ or faster clock speed.

[0039] The operator can use interface device 219 to interact with thesystem 200, such as to indicate which portion of the scanned body isdesired to be explored. The interface device 219 can further be used tocontrol the image being displayed, including the angle, size, rotation,navigational position and the like.

[0040] Scanning device 205 and terminal 216, or parts thereof, can bepart of the same unit. Numerous CT, MM, and US systems are suitable forsuch applications. A single platform may be used to receive the scannedimage data, convert the image data to 3D voxels if necessary and performthe guided navigation.

[0041] In many virtual visualization and surgical planning operations,the use and registration of multiple images is desirable. This can bethe result of acquiring multiple sets of image data from a region overtime, acquiring multiple image sets of a region using different patientorientations, or by acquiring multiple image sets of a region usingvarious technologies to exploit the advantages of various imagingtechniques. When multiple sets of image data are acquired, such imagesets must by registered in some manner in order to reasonably use thedata in each image set. Exact image registration, for example, a voxelspecific registration process, can be very difficult to achieve and isprocessor intense. In many applications, such precise registration isnot required.

[0042]FIG. 3 is a flow diagram illustrating an overview of a flexibleregistration process for aligning multiple sets of image data, such asin a series of single-modality images and in multi-modality imagingwhich can be practiced using the system of FIG. 2 and variationsthereof. This process will be described in connection with thevisualization and analysis of plaque build up within a blood vessel, butis extendible into many forms of single-modality image sets,multi-modality imaging and other multiple image set registrationapplications. In step 305, a first set of image data is acquired. In thecase of arterial stenosis, the first set of image data may be acquiredusing 3D time-of-flight magnetic resonance angiograpy (TOF MRA). The 3DTOF MRA has proven to be an effective imaging technique for providinggood enhancement of the arterial lumen. Thus, the gross structure of thearterial lumen and the suspect location of a plaque build up can beidentified using this imaging technique. However, TOF MRA does notprovide optimal resolution at distinguishing the characteristics amongsoft tissues and the components of the plaque.

[0043] In step 310, a second set of image data is acquired. As notedabove, this can be image data acquired at a different time, a differentpatient orientation or with a different imaging technology. In the caseof arterial stenosis, a second imaging mode is used to provide a set ofimage data which can be segmented to distinguish characteristics of thesoft tissue and plaque. For example, T₁-weighted MR imaging (T₁-MRI) hasbeen found to provide good contrast and can be performed in the same MRIseries as the TOF MRA. In order to obtain additional information, yetanother set of image data can be acquired to provide features not easilydistinguished in the previous image data sets. For example, MRI datadoes not readily distinguish a calcification component of the plaque. Toproperly determine the total volume and composition of the plaque, a CTimage dataset of the region can also be acquired.

[0044] Using one set of the image data acquired in steps 305 and 310 themajor structural characteristics of the region, including possiblesuspect regions, can be ascertained (step 315). As noted above, the TOFMRA image data provides a good rendering of the 3D arterial lumen andcan be used to identify any regions of possible stenosis.

[0045]FIG. 4 is a pictorial diagram illustrating a section of a bloodvessel with plaque build up resulting in stenosis. The arterial lumen400 is shown with a region of plaque 405 impinging on the lumenresulting in stenosis. The lumen can be characterized by a diameter ofthe normal vessel (B) 410 and a minimum diameter at the location of thestenosis (A) 415. The degree of stenosis can be characterized by thedefinition set forth by the North American Symptomatic CarotidEndarterectomy Trial (NASCET). The NASCET is determined as:$\frac{B - A}{B} \times \quad 100\quad \%$

[0046]FIG. 5 is a pictorial diagram illustrating a branched section of ablood vessel with plaque build up resulting in stenosis. From step 315,the lumen of the artery 505 can be extracted and the regions of possiblestenosis 505 can be identified. In addition, a centerline skeleton ofthe arterial structure can also be determined. In FIG. 5, the arteryexhibits a first arterial section 510 which branches at a bifurcationpoint 515 into an internal arterial branch 520 and an external arterialbranch 525.

[0047] After at least one set of image data has been processed todetermine the structural characteristics of the region (step 315), atleast one control point is established using a consistently identifiablestructural characteristic of the region (step 320). Referring to FIG. 5,in the case of a branched arterial structure, the bifurcation point 515can be determined and selected as the control point for subsequentregistration operations.

[0048]FIG. 6 is a pictorial diagram of a branched section of a bloodvessel further illustrating the location of the bifurcation point whichcan be used as a control point for image data registration. From theskeletal centerline of the arterial lumen the bifurcation point can bedetermined. Referring to FIG. 6, the bifurcation point is refined asbeing the center of a maximal diameter sphere that fits within the lumenat the bifurcation region. Using the center of the maximal sphere as theorigin, a coordinate system, such as X, Y and Z axes, can be defined inthe region for subsequent registration. A number of different coordinatesystems can be used. For example, a Cartesian coordinate system can beused which corresponds to the orientation of the human body beingimaged, e.g., with the Z-axis aligned along the height of the body(e.g., from head to toe), the Y-axis oriented from back to front and theX axis running laterally (e.g., from left to right). The units of lengthfor the coordinate system can be arbitrarily set to one voxel or someother convenient unit of measure. The absolute magnitude of the unitsused will vary based on the acquisition properties of the imagingscanner being used.

[0049]FIG. 7A further illustrates the determination of the control pointand the registration between a spherical model and its deformable model.As illustrated in FIG. 7A, three mutually orthogonal planes alignedalong a Cartesian coordinate system are defined which can readily beapplied to a deformable model, such as that illustrated in FIG. 7B. Thetwo coordinate systems can be registered by aligning them each to areference coordinate system. The reference coordinate system can be onewhich is associated with anatomical axes of the human body. Thisflexibility allows registration of multiple image sets despitevariations in resolution and other variables.

[0050] Returning to FIGS. 3 and 5, using the control point 515 and thestructural characteristics of the region, the approximate location ofsuspect regions is determined in at least one of the sets of image data(step 325). In the case of carotid stenosis, the location of the suspectregions is determined by finding the areas of minimum diameter 530 inthe lumen, which correspond to areas of maximum stenosis in a region ofplaque. The location of the suspect region can be determined in a numberof ways. For example, the curve length L from the control point 515 tothe point of minimum diameter 530 along the skeletal centerline providesa first measure of location. In addition, the “vertical” distance H,which is measured along the imaginary Z-axis of the coordinate systemframing the control point 515, can be used to determine the relativeposition of the suspect region with respect to the control point.

[0051] After the control point(s) and the position of suspect areas havebeen identified in a first image data set, registration of the firstimage data set with the additional image data sets can take place (step330). The process of flexible registration entails identifying thecontrol point(s) in each image data set and using the control points asa datum for registration of the image data sets. For example, in theexample of a branched arterial structure discussed above, thebifurcation point would be located in each set of image data. In thecase of T₁-weighted MR images, generally the bifurcation point ismanually identified by a user of the system, such as a physician. Theimage data sets are registered to one another using the control point(s)and are aligned to the coordinate system which is assigned. Assuming anaffine model of the region being imaged, three control points can beused to fully automatically register the coordinate systems of the imagedata sets in 3D space. With less than three control points, theregistration process can be initially automated and then completed usingmanual intervention to optimize the registration process.

[0052] Using an appropriate transfer function, the voxels of thearterial lumen can be displayed as a translucent structure with theskeletal centerline therein. In addition, the suspect regions can behighlighted on the skeletal centerline by emphasizing the suspectregions, such as by using a different color to display suspect regions.

[0053] To determine the location of stenosis, the information from thefirst image data set, in this case the previously processed results ofTOF MRA, are preferably used. The bifurcation point is detected in theT₁-weighted MR images, generally, manually by physician. The samecoordinate system defined with respect to the first image data set isthen applied to the T₁-weighted MR images. Starting at the controlpoint, the vertical distance H is measured along the Z-axis to arrive atthe position corresponding to the approximate location of the stenosisin the second set of image data (step 335). In the T₁-weighted MR imagesthe vessel lumen is difficult to differentiate. However, if an imagingtechnology is employed which offers adequate delineation of the vessellumen, the distance L along the center line of the lumen to the suspectregion can also be used. Exact voxel accuracy of location is notnecessary in aligning the images. The important criteria is that asuspect region in the image data covers the entire region of stenosisand plaque which has been identified in step 325.

[0054] Image segmentation can then be performed in the region of plaqueto determine the nature and composition of the plaque build up (step340). If the plaque has fatty components, those components will exhibita higher intensity in the T₁-weighted MRI image as compared to both thevessel lumen and the vessel wall. Such components can be detected andidentified in the segmentation results of the T₁-weighted MR images. Thevolume of the fatty components can be quantified and a 3D model of theplaque can be created.

[0055] As noted above, the TOF MRA data set provides good delineation ofthe lumen whereas the T₁-weighted MR data provides good segmentationresults for the fatty components of the plaque causing the stenosis. Byemploying flexible image registration and then fusing the registereddatasets into a single display data set, a single display providing theenhanced features of both the arterial structure and the plaque can beprovided (step 345).

[0056] Image fusion can take on several forms but generally results inat least a portion of the features of each image set being presentedsimultaneously or at least selectively on a common display and in acommon registration. Generally, each set of image data, whethermulti-modality image data or a series of single-modality image datasets, is pre-processed through an image segmentation process. Imagesegmentation is used to emphasize and extract various features in thesets of image data. One or more sets of image data can be selected todisplay the raw image data. In the case of arterial imaging for carotidstenosis, this would generally by the TOF MRA image data or US imagedata which effectively delineates the arterial lumen. Segmented portionsof the other image data sets can be simultaneously or selectivelyrendered on the display. Each of the sets of image data can be renderedas a layer in the display to effectively render the region of interest.For example, in viewing a region of stenosis, an operator canselectively heighten the image data from the T₁-weighted MR data toemphasize fatty components of the plaque or the CT image data toemphasize the calcification components of the plaque. The rendering canbe performed initially in an automated manner than optimized by theuser, generally by use of a GUI interface. In addition, the operator cangenerally navigate in the image data, change views, zoom, etc. Theprocessing of the present method effects these navigation andmanipulation commands on all fused image data sets such that allinteraction is effected on all image models. This is possible becausethe sets of image data are registered to one or more common controlpoints and are aligned to a common coordinate system.

[0057] The T₁-weighted MR data provides good segmentation results forthe fatty components of the plaque but generally can not reliablydistinguish the calcification components of the plaque. Generally, CTimage data exhibits higher spatial resolution and better contrastbetween calcification and other tissues. Thus, it is desirable to scanthe region with a CT scanner to acquire CT image data which can also beflexibly registered with the TOF MRA data set.

[0058] In the CT image data set, the vessel lumen and wall have similarintensity values as that of soft tissue. As a result, the arterial lumenis difficult to delineate in a CT image. A first approach to overcomethis issue is to perform CT image scanning twice. In a first imagingscan, CT angiography (CTA) is acquired following the injection of acontrast agent in the blood. From the CTA image data, the carotid lumenwill be delineated. The skeleton along the centerline of the lumen canbe generated and the bifurcation point can then be determined. Bytracing the skeleton up-stream the distance L, the stenosis location canbe detected. The second CT imaging scan is acquired with the samescanning settings, but without the presence of the contrast agent. Inthe second scan, only calcification regions are depicted as regions ofincreased intensity. Since the field of view (FOV) and scanning settingare identical, the two scans are assumed to be registered atvoxel-specific accuracy. The region of stenosis in the second scan isthen considered to be the same as that in the CTA scan. Then, a suitablesegmentation algorithm can be applied to the region of stenosis toextract the region of calcification. The volume of the calcification canbe measured and a 3D model of it can be created.

[0059] The second approach to using the CT imaging data to determine theregions of calcification does not use a CTA scan with a contrast agent.Instead, the stenosis region is detected manually in the CT images. Thistask is generally performed by a physician with help of image renderingtool. Various forms of image rendering tools are known in the art andgenerally provide a graphical user interface (GUI) which allows the userto readily manipulate both the individual image slices and the 3Dvirtual model on a computer screen. An example of one such GUI displayis depicted in FIG. 10. Anatomical knowledge of the structure, such asthe carotid artery, is generally used in this procedure.

[0060]FIG. 10 illustrates an exemplary GUI display of a fused display ofa segmented object 1005 within the raw CT image data of the region ofinterest 1000, in this case the bladder of a patient. In the viewdepicted in FIG. 10, the image data is presented as a 3D body in 2Dforms, i.e. the image is shown in cut away slices in two orthogonaldirections and projecting in a third direction. The GUI includeson-screen controls 1010 which can be selected and modified using eithera computer keyboard or digital pointer, such as a mouse, trackball andthe like. Using the GUI the user can inspect details of the image, editthe segmentation results, change the view and the like.

[0061] After detecting a suspect region, such as an area of possiblestenosis, the suspect region image is segmented and processedsubstantially as described above. While this second process requiresfurther manual intervention, it does not require the use of a second CTscan or the use of a contrast agent.

[0062] While in the above example, TOF-MRA and CTA imaging weredescribed as imaging methods to determine the arterial lumen, it will beappreciated that US imaging can also be used in place of one or both ofthese imaging techniques.

Determination of Likelihood of Rupture

[0063] The plaque which forms in the arterial lumen and results instenosis is generally composed of blood volume, fatty components andcalcification components. The total volume of the plaque can be obtainedby summation of the volume of fatty components identified in theT₁-weighted MR image data whereas the calcification components can bedetermined from the segmented CT image data. The percentage of eachcomponent can be calculated. The larger the volume of blood or itspercentage, the higher the probability of a rupture of the plaque, withincreased attendant risk of stroke. In addition, a close proximity ofblood deposits near the surface of the plaque build up also provides anindication of high risk of rupture. The level of risk is notquantifiable in absolute terms. However, by employing an arbitrary scalefor both blood volume and surface depth, a weighted relative risk factorcan be associated with the stenotic region. Based on such “risk scales,”regions of risk can be highlighted, such as by use of one or moredifferent display colors, to alert the user of the risk. The risk factorvalues can also be compared to predetermined or user defined thresholdvalues, which if met or exceeded, would trigger more active alarmindications to the user.

[0064] The volume of the blood inside the plaque and its location nearthe boundary are two indicators for the likelihood of rupture.Therefore, accurate segmentation of the plaque components andidentification of their spatial relation are objectives in stroke riskassessment. In the present systems and methods information from (1)CT/CTA (calcification), (2) TOF MRA or US(lumen), and (3) T₁-weightedMRI (fatty, blood, and others) images are simultaneously used to renderthis information. The information is represented together to show thetissue spatial relation in 3D space. This requires the use of both imageregistration and image segmentation computer processing techniques. Eachof these operations benefit from a degree of user editing. The editingcan be performed by displaying the segmented objects inside the rawimage or by displaying all segmented objects together, as depicted inFIG. 10.

[0065]FIG. 8 is a flow chart summarizing the application of the methodof FIG. 3 as specifically applied to the visualization and analysis ofan arterial structure for carotid stenosis. In step 805 image data isacquired to visualize the arterial wall and lumen. This generally takesthe form of TOF MRA or US image data. After image segmentation, thearterial wall and lumen are clearly distinguishable in the MRA or USimage data (step 810). From the MRA image which illustrates the lumen,the centerlines of the lumen are calculated and the bifurcation point islocated as the control point (step 815). Suspect regions are located inthe lumen which may reflect areas of stenosis (step 820). The positionof these suspect regions is then determined with respect to the controlpoint and an arbitrary coordinate system which is defined around thecontrol point.

[0066] In addition to the TOF MRA image data, a second set of image datais acquired to highlight fatty components within the suspect regions(step 825). Generally, a T₁-weighted MRI image data set has been foundacceptable for this purpose. The control point and coordinate systemwhich was defined in the first image data are then located and placed inthe second image data set and the position of the suspect regions arelocated (step 830). Image segmentation is performed on the second imagedata set in the suspect regions to highlight the fatty components of theplaque deposits.

[0067] A third image data set is acquired to delineate the calcificationcomponents of the plaque deposits (step 840). This third image data setcan take the form of CT image data. Within the third image data set, thecontrol point and coordinate system are identified and placed and theposition of the suspect regions are identified (step 845). Again, imagesegmentation in the suspect regions can then be performed to delineatethe regions of calcification (step 850).

[0068] Using the control point and coordinate system which has beenlocated in each set of image data, the sets of image data are registeredwith one another and fused into a single display image which can bepresented on the display of the visualization system (step 855). Usingthe fused image, the user can readily analyze the suspect regions, forexample to determine the volume and composition of any plaque deposits,the likelihood of rupture and the quality of the arterial surface forreceiving a stent graft (step 860).

Spatial Inhomogeneity Correction

[0069] In MR image data, the inter-slice and intra-slice image dataoften contains spatial inhomogeneities which adversely affect automatedanalysis of these images. Correction for spatial intensityinhomogeneities in inter- and intra-slices provides improved results inautomated quantitative analysis of MR images. The inhomogeneitieslocally alter image intensity mean, median, and variance. It is knownthat they vary over time and with different acquisition parameters andsubjects. For virtual treatment planning, a general purpose correctionalgorithm that can work for any scanning protocol and anatomy region isdesirable.

[0070] A correction method based on known renormalization group theoryhas now been found to be useful for the spatial inhomogeneity of MRimage data. As illustrated in FIG. 9, following image data acquisition(step 905), the correction method consists of two primary steps. Thefirst step is to estimate an inhomogenous bias field (step 910). Thesecond step is to correct the MR image based on that bias field (step915).

[0071] The notations which will be used in connection with the presentrenormalization transformation (RT) method for MR image are nowpresented. Let S be a single-connected, finite subset of the 2D integerlattice. One can consider S as the pixel locations of a region ofinterest (ROI) in a 2D image. Let X={x_(i): iεS} be an image, wherex_(i) is the intensity value of pixel i and it takes value within afinite set. Denote ∂_(i) and |∂_(i)| as the first order neighborhood ofthe pixel i and the number of pixels inside this neighborhood (pixel iitself is also a member of this neighbor). The RT is defined as:${T(X)} = {\left\{ {y_{i}:{{y_{i} \equiv {\frac{\sum\limits_{k \in {\partial_{i}{\bigcap S}}}^{\quad}\quad x_{k}}{\left| {\partial_{i}{\bigcap S}} \right|},\quad i}} \in S}} \right\}.\quad 2.1}$

[0072] T(X) is a new image with the same pixel location and itsintensity value of each voxel is the average intensity of those ofpixels in its neighborhood. The RT is a linear transformation and can becalculated locally. This makes it possible for a fast parallelimplementation.

[0073] The estimation of the bias field is to apply RT iteratively tothe MR image X. When the condition$\max\limits_{i \in S}\left| {{T^{(n)}(X)}_{i} - {T^{({n + 1})}(X)}_{i}} \middle| {\leq t} \right.$

[0074] is satisfied, T^((n))(X) is the estimation of the bias field,where n is the number of iteration and t is a pre-set threshold. Forexample, we can set t as the average intensity of the image X dividingby 3000.

[0075] The bias field is often modeled as the following classical way:

Y=X·f+n  2.2

[0076] where Y is the obtained MR image, X is the ideal image withoutinhomogeneity, f is the bias field, and n is the noise. Ignoring thenoise component, the image data is the multiplication of the ideal imageand the bias field. Following this model, the correction procedure canbe performed by dividing the image data by the estimated inhomogeneitybias field. However, in the RT model, the RT is a linear transformationrather than a multiplication transformation. Therefore, the followingmodel is more accurate for RT method:

Y=X+f+n  2.3

[0077] Where Y, X, n, and f were defined above. For correctionalgorithms ignoring the noise, we induce the following model:

X=Y+α·({circumflex over (f)}−M ₀)  2.4

[0078] where α({circumflex over (f)}−M₀) is a linear transform appliedon the estimated bias field {circumflex over (f)}, α is a scaling factorand M₀ is an image calculated from X. This model assumes that all pixelsrepresent a common tissue type. In practical situations, however, thisis not the case. Hence, a correction model is applied in a flexible wayto ensure that the assumption is satisfied to the greatest degreepossible. For example, the image can be processed on a row by row orcolumn by column basis. Alternatively, image segmentation can beperformed initially and regions of substantially common tissue types canbe delineated and separately processed by applying the model to eachregion. The determination of α and M₀ should also be adaptive to thepractical situation. For example, we can set α=−1 and compute M₀ as theaverage intensity of the image.

Image Segmentation

[0079] A difficulty encountered in imaging is that several of therelevant anatomical structures have similar intensity values on the CTimage. This can make it difficult to distinguish the various structureswhich overlap or are interrelated. To address this problem, a two-levelimage segmentation process can be employed. The two-level segmentationprocess involves low-level processing of the voxels in the region ofinterest followed by high-level organ extraction. During the low-levelprocessing, the voxels of the 3D dataset are clustered into groups basedon an intensity feature of the voxels, which can be measured by anassociated local intensity value vector. This can be determined using amodified self-adaptive on-line vector quantization algorithm, such as isdescribed in the article “A self-adaptive on-line vector quantizationalgorithm for MRI segmentation,” by Chen et al. in the proceedings ofThe 7th Scientific Meeting of ISMRM, May 1999, Philadelphia, which ishereby incorporated by reference. In the low-level classification, eachvoxel is associated with a local vector which is defined in 3D space.From the local vectors, a feature vector series can be derived using acomponents analysis which is well known in the art. The feature vectorsare then clustered using a self-adaptive on-line vector quantizationalgorithm. The voxels are then grouped according to the classificationof their feature vectors and are assigned an integer value representingthis classification.

[0080] After the low-level processing is complete, the high level organextraction processing can follow. Initially, a user locates a seed, orstarting point, within regions representative of soft tissue, bone andair spaces. The system then applies a region growing algorithm startingfrom the seed points to extract the anatomical features of the region.

[0081] Bone structure, which presents different contrast compared to thesurrounding tissues is fairly easy to automatically segment. However,certain structures may require additional user input to fully delineatethese structures. For example, the soft tissue of the inner ear presentsa similar intensity value on CT images as compared to the surroundingsoft tissue. Thus, to insure proper extraction of such features, it maybe desirable for the user to manually delineate the outline of thestructure by manually tracing the contour on one or more of the imageslices.

[0082] While the above described two level image segmentation ispreferred, any method which provides accurate delineation of theneighboring structures in a region of interest can be used in thepractice of the present treatment planning method. One such technique isdescribed in the article “On segmentation of colon lumen for virtualcolonoscopy” by Liang et al., Proceedings of SPIE Medical Imaging, pp270-278, February 1999, San Diego.

[0083] Once image segmentation is performed, 3D image generation can beperformed for each of the segmented objects using a number of knowntechniques, such as the Marching Cubes algorithm, which reconstructs theouter polygonal surface. However, because of the complexity of manyanatomical structures, interactive rendering of all polygons in thedisplay for each change to a portion of the display is processor intenseand unduly costly. As more colors and surfaces are delineated in thedisplayed image, this burden increases. To minimize the processingoverhead, the volume image dataset can be stored in a partitioned datastructure, such as a binary space-partitioning (BSP) tree, in which thelarge dataset is parsed into relatively small portions which are storedin leaf nodes of the data structure. By identifying which leaf nodes areeffected by any given operation, and only performing operations, such asConstructive Solid Geometry (CSG) operations, on the effected leafnodes, the processing burden for interactive operations can besignificantly reduced. As will be set forth in more detail below, theprocessing burden can be further reduced by use of a level of detail(LOD) rendering mode and/or a wavelet transformation to reduce the datavolume.

[0084] As noted above, when a large dataset is involved, it may berequired, or at least desirable, to reduce the size of the dataset tospeed up processing and reduce processing cost. Noting that the treestructure can be preserved within a range of scales, the large volumecan be shrunk to a smaller scale space for structure analysis.

[0085] A shrinking method based on multiresolution analysis theory canbe used. The input data is the stack of binary images of the same sizewhich can be obtained from the segmentation results of the CT or MRIscan. The x-direction is taken along the slice image width, they-direction is along the slice image height, and the z-direction isalong the direction of slice by slice. The foreground voxels in the treevolume are set to value of 128 (maximum) and the background voxels areset to value 0 (minimum). A Daubechies' bi-orthogonal wavelet transformwith all rational coefficients can be employed. This one-dimensional(1D) discrete wavelet transformation (DWT) is first applied along to thex-direction row by row. From application of the DWT only the lowerfrequency components are retained and packed. The computation ispreferably implemented in floating points. Noting that the DWT isapplied to the binary signal, there are two kinds of nonzerocoefficients which result in the lower frequency component. The first isof value 128 and this kind of coefficient is located in the interior ofthe volume. The second is of a value not equal to 128 and thesecoefficients locate the boundary of the volume.

[0086] The coefficients of the second kind are compared against apredetermined threshold value. If the absolute value of the coefficientsis larger than a pre-set threshold T1, the value of the coefficient isset to 128; otherwise, it is set to 0. This results in a stack of binaryimages with a row size of half of the original dataset. The same DWT isthen applied to the resulting dataset along the y-direction column bycolumn, where the similar thresholding is employed to the lowerfrequency components. The result is again a stack of binary images, butnow with both half row and column size as compared to the originaldataset. Finally, the DWT is applied to the last result along thez-direction and the lower frequency components are retained. This stepcompletes the first level decomposition.

[0087] The resulting dataset of the first level decomposition is of halfsize in all three directions as compared to the original dataset. If theshrinking procedure stops at this level, the finial thresholding isapplied. It revalues those coefficients of nonzero and non-128 value. Ifthe absolute value of this kind of coefficient is larger than a pre-setthreshold T2, it will be revalued as 128; otherwise, it is revalued as0. If further shrinking is needed, the same thresholding algorithm isapplied with the threshold T1. Further shrinking proceeds as previouslydescribed, but is applied to the dataset shrunk at the last previouslevel. The decomposition procedure can be recursively applied until theresulting volume meets the desired reduced data volume. In the casewhere the slice images are of 512×512 pixel size, the maximumdecomposition level is usually three, resulting in a 64×64 reduced pixelsize.

[0088] The volume is isotropically shrank in all directions with thepresented method. The two pre-set thresholds, T1 and T2, are used tocontrol the degree of shrinking. If the volume is significantly overshrunk, connectivity may be lost in the reduced volume. If it is overshrunk to a lesser degree, two separate branches may merge into onebranch in the reduced volume dataset. The larger the two thresholdvalues, the thinner the reduced volume is. The range of those twothresholds is [0, r×128], where 0<r<1. Preferably, the range for virtualendoscopy is rε(0.08, 0.28) for T1 and rε(0.7, 0.98) for T2. The exactdetermination is dependant on the feature size of the particularapplication and is selected to achieve reduction while retaining thefidelity of the structure information in the shrunk volume.

[0089] The flexible image registration techniques and spatialinhomogeneity corrections are applicable to surgical planning processes.For example, in the case of angioplasty, the surgeon inserts a catheterinto an occluded artery and inflates a balloon at the end of thecatheter to force the occluded artery open and to expand a stent whichmaintains the opening. While this has become a common procedure, it isnot without risk. For example, the arterial occlusion is generallyrelated to a build up of plaque and fatty deposits on the arterial wall.If a portion of these deposits are dislodged during the angioplastyprocess or if the region of plaque ruptures, there is a risk of strokeand other complications. Using the present method of treatment planning,the artery can be imaged and, through image segmentation, the quantityand nature of the plaque deposits can be determined. The severity of theocclusion can be viewed by the surgeon who can navigate in the 3D imagewithin the artery. A virtual intervention can then be performed, i.e.,placing a virtual catheter within the arterial volume and expanding avirtual stent, and the results observed. If problems are observed, suchas a high risk of rupture, the user can then alter the course oftreatment to minimize the risk. The virtual catheter and stent require adynamic model that conforms to the contours of the interior surface ofthe arterial wall. Such a model is analogous to the force field modelpreviously used in guiding a virtual camera along a fly path inperforming virtual colonoscopy.

[0090] Known volume rendering techniques use one or more definedtransfer functions to map different ranges of sample values of theoriginal volume data to different colors, opacities and otherdisplayable parameters for navigation and viewing. During navigation,the selected transfer function generally assigns maximum opacity to thewall of the object being viewed. However, once a suspicious area isdetected during virtual examination, a user, such as a physician, caninteractively change the transfer function assigned during the volumerendering procedure such that the outer surface being viewed becomessubstantially transparent, allowing the interior structure of the regionto be viewed. Using a number of predetermined transfer functions, thesuspicious area can be viewed at a number of different depths, withvarying degrees of opacity assigned throughout the process. In addition,the shape of the region and texture of the region undergoing virtualexamination can be analyzed to determine the nature of the region, suchas a likelihood of cancerous tissue residing in a region being biopsied.

[0091] Although the present invention has been described in connectionwith specific exemplary embodiments, it should be understood thatvarious changes, substitutions and alterations can be made to thedisclosed embodiments without departing from the spirit and scope of theinvention as set forth in the appended claims.

1. A computer based method of visualizing a region using multiple imagedata sets comprising: acquiring first volumetric image data of a region;acquiring at least second volumetric image data of the region;determining structural features of the region using the first volumetricimage data determining at least one control point in the region using atleast one structural feature in the first volumetric image data;identifying the at least one control point in the region in the at leastsecond image data of the region; registering the first image data andthe at least second image data using said at least one control point;and fusing the registered first image data and at least second imagedata into a common display data set.
 2. The method of visualizing aregion using multiple image data sets of claim 1, wherein the regionincludes a branched structure having a common section and at least twobranches extending from the common section forming a bifurcation regionand wherein said control point is the bifurcation point of the branchedstructure.
 3. The method of visualizing a region using multiple imagedata sets of claim 2, wherein the bifurcation point is determined as thecenter of an imaginary maximal diameter sphere which is fit within thebifurcation region.
 4. The method of visualizing a region using multipleimage data sets of claim 3, wherein the branched structure in a bloodvessel.
 5. The method of visualizing a region using multiple image datasets of claim 1, wherein the at least second volumetric image data isacquired using a mode which is different from a mode used to acquire thefirst volumetric image data.
 6. The method of visualizing a region usingmultiple image data sets of claim 5, wherein the mode used to acquirethe at least second image data is selected to provide complimentaryimaging characteristics to the mode used to acquire the first imagedata.
 7. The method of visualizing a region using multiple image datasets of claim 1 wherein the second volumetric image data is acquired ata substantially different time from said first volumetric image data. 8.The method of visualizing a region using multiple image data sets ofclaim 1, wherein the second volumetric image data is acquired with thepatient in an orientation which is different than an orientation usedduring acquisition of the first volumetric image data.
 9. A method ofvisualizing and analyzing a region of artery comprising: acquiring firstvolumetric image data of a region of artery using a first acquisitionmode; acquiring at least second volumetric image data of the region ofartery using a second acquisition mode, said second acquisition modebeing different than the first acquisition mode; determining structuralfeatures of the artery using the first volumetric image data;determining a bifurcation point in the artery using the first volumetricimage data; defining a coordinate system with the bifurcation point asthe origin; identifying the bifurcation point in the artery in the atleast second image data of the region and applying the definedcoordinate system to the at least second volumetric image data;registering the first volumetric image data and the at least secondvolumetric image data using the bifurcation point and coordinate system;and fusing the registered first volumetric image data and at leastvolumetric second image data into a common display data set.
 10. Themethod of visualizing and analyzing a region of artery of claim 9,wherein the structural features include at least one region of stenosisresulting at least in part from a plaque deposit, and wherein the methodfurther comprises the steps: identifying the location of the at leastone region of stenosis in the first volumetric image data; using thelocation from the first volumetric image data, determining the locationof the stenosis in the at least second volumetric image data.
 11. Themethod of visualizing and analyzing a region of artery of claim 9,wherein the step of identifying the location of the at least one regionof stenosis includes measuring the distance along at least one axis ofthe coordinate system from the bifurcation point to a point in theregion of stenosis
 12. The method of visualizing and analyzing a regionof artery of claim 10, wherein the at least second acquisition mode isselected to delineate fatty components of the plaque deposits.
 13. Themethod of visualizing and analyzing a region of artery of claim 10,wherein the at least second acquisition mode is selected to delineatecalcification components of the plaque deposits.
 14. The method ofvisualizing and analyzing a region of artery of claim 10, wherein the atleast second acquisition mode includes second and third acquisitionmodes, at least one of the second and third acquisition modes beingselected to delineate fatty components of plaque deposits and at leastone of the second and third acquisition modes being selected todelineate calcification components of plaque deposits.
 15. The method ofvisualizing and analyzing a region of artery of claim 14, wherein thefirst acquisition mode is Magnetic Resonance Angiograpy.
 16. The methodof visualizing and analyzing a region of artery of claim 14, wherein thefirst acquisition mode is ultrasound imaging.
 17. The method ofvisualizing and analyzing a region of artery of claim 14, wherein thesecond acquisition mode is T₁-weighted MR imaging.
 18. The method ofvisualizing and analyzing a region of artery of claim 14, wherein thesecond acquisition mode is T₂-weighted MR imaging.
 19. The method ofvisualizing and analyzing a region of artery of claim 14, wherein thefirst acquisition mode is time-of-flight Magnetic Resonance Angiograpy,the second acquisition mode is T1-weighted MRI imaging, and the thirdacquisition mode is Computed Tomography imaging.
 20. The method ofvisualizing and analyzing a region of artery of claim 9, furthercomprising analyzing a plaque deposit identified in the fused imagedata.
 21. The method of visualizing and analyzing a region of artery ofclaim 20, wherein the analyzing step includes determining thecomposition of the plaque deposit.
 22. The method of visualizing andanalyzing a region of artery of claim 21, wherein the analyzing stepincludes determining the volume and position of fatty components in theplaque deposit.
 23. The method of visualizing and analyzing a region ofartery of claim 21, wherein the analyzing step includes determining thevolume and position of calcification components in the plaque deposit.24. The method of visualizing and analyzing a region of artery of claim21, wherein the analyzing step includes determining the volume andposition of blood components in the plaque deposit.
 25. The method ofvisualizing and analyzing a region of artery of claim 20, wherein theanalyzing step further comprises assessing the risk of rupture of theplaque deposit.
 26. A method for correcting spatial inhomogeneity inmagnetic resonance imaging data comprising: acquiring MRI image data;estimating a bias field associated with the MRI image data; and applyingthe bias field to the MRI image data to provide corrected MRI imagedata.